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Abstract 



We consider the problem of pricing path-dependent contingent claims. Classically, this problem 
can be cast into the Black-Scholes valuation framework through inclusion of the path-dependent 
variables into the state space. This leads to solving a degenerate advection-diffusion Partial 
Differential Equation (PDE). Standard Finite Difference (FD) methods are known to be inade- 
quate for solving such degenerate PDE. Hence, path-dependent European claims are typically 
priced through Monte-Carlo simulation. To date, there is no numerical method for pricing 
path-dependent American claims. 

We first establish necessary and sufficient conditions amenable to a Lie algebraic characteriza- 
tion, under which degenerate diffusions can be reduced to lower-dimensional non-degenerate 
diffusions on a sub-manifold of the underlying asset space. We apply these results to path- 
dependent options. Then, we describe a new numerical technique, called Forward Shooting 
Grid (FSG) method, that efficiently copes with degenerate diffusion PDE. Finally, we show 
that the FSG method is unconditionally stable and convergent. 

The FSG method has been implemented for a number of popular path-dependent options, and 
proved to be much faster than traditional Monte Carlo simulation, for a comparable accuracy. 
Depending on the type of option, the computation time lies between 1 and 1 5 seconds on a PC, 
for a 0.1 % precision. 

The FSG method is also the first capable of dealing with the early exercise condition of American 
options. Furthermore, when the stock price S follows a binomial process, the method computes 
the exact price of any American lookback option on 5. The same is true for barrier options, 
such as up-and-in or down-and-out options. 

Several numerical examples are presented and discussed, showing in particular that the Snell 
envelope upper bounds obtained on American lookback prices are quite overestimated, and 
also that the usual geometric average approximation to the arithmetic average-rate option is 
fairly inaccurate. 
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1 Introduction 

Path-dependent options are options whose payoffs depend on historical values of the underlying 
asset over a given time period as well as its current value. Well-known examples are lookback 
call (put) options, which give their owners the right to buy (sell) the asset at an exercise price 
equal to the minimum (maximum) price of the asset over the life of the option. Many other 
variants exist, e.g. capped options, barrier options etc. 

Average-rate or Asian options constitute another class of path-dependent instruments, whose 
payoffs depend on the arithmetic average value of the asset price for some time period. Popular 
examples include fixed-strike, and floating-strike average-rate options. 

Since their introduction in 1982, path-dependent options found their way in several places 
such as common stocks and foreign exchange markets, by meeting specific risk management 
and investment needs. As evidence of their increasing popularity, some instruments such as 
capped options have recently begun to be traded on the Chicago Board of Options Exchange 
and American Stock Exchange, see e.g. Hunter and Stowe (1992). 

The modern approach to path-dependent option pricing relies on the dynamic hedging principle 
of the Black-Scholes model (Black and Scholes (1973)). In their seminal work, Goldman et al. 
(1979) have shown that a hedge portfolio could be constructed for an option to buy at the 
historical maximum, and that closed-form valuation formulas exist in the European case. 
Garman (1989) developed another valuation model, which separates the lookback option into 
two underlying options, and gives furthermore the ability to price a European option on an asset 
that pays dividends. The Asian option is analyzed by Yor (1989). Conze and Viswanathan 
(1991) derive explicit valuation formulas of most European barrier options, as well as some 
upper bounds in the American case. Kemna and Vorst (1990) use Monte Carlo simulation with 
a specific variance reduction method to compute the price of fixed-strike average-rate options. 

The path-dependent pricing problem can be cast into the classical Black-Scholes valuation 
framework through inclusion of the path dependent variables into the state space (see e.g. 
Stanton (1989)). In a few simple cases, the resulting augmented PDE admits a closed form 
solution. However, the use of numerical techniques is mandatory in most practical situations. 
The augmented PDE associated with a path-dependent problem is generally degenerate, i.e. 
the instantaneous covariance matrix is singular. Finite Difference (FD) methods are numer- 
ically stable, but they introduce a spurious additional numerical diffusion, and therefore do 
not converge towards the theoretical solution. The only practical technique to date consists in 
computing the price through Monte Carlo simulation by means of the Feynman-Kac represen- 
tation. This approach is satisfactory for pricing European path-dependent contingent claims. 
However, it cannot deal with the early exercise conditions of American claims. 

In this paper, we first establish necessary and sufficient conditions under which degenerate 
diffusion PDE can be reduced to lower-dimensional non-degenerate PDE on a sub-manifold of 
the state space. A degenerate PDE that can be reduced in such a way is called holonomic. It is 
called non-holonomic otherwise. It is important to determine whether a given PDE is holonomic 
or not, since this property influences both the type of numerical integration technique that can 
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be used, and the memory required to conduct the integration. 

For this purpose, we use the concept of symmetric multiplication for stochastic integrals, 
which leads to the notion of Fisk-Stratonovitch differential. Then, we use standard results 
regarding discrete approximations of multidimensional diffusion processes that establish a link 
between stochastic differential equations and deterministic optimal control theory. Finally, we 
use standard results on the accessibility of deterministic control systems. This leads to a Lie 
algebraic characterization of holonomic PDE. 

Then, we apply these results to popular types of path-dependent pricing problems. We show 
in particular that the average-rate option pricing PDE is non-holonomic. Finally, we present a 
new numerical technique, called Forward Shooting Grid method, that efficiently copes with the 
degeneracy of these non-holonomic PDE. In particular, we show that, unlike FD methods, the 
FSG method is unconditionally stable and convergent in the presence of arbitrary degeneracies. 

The FSG method has been implemented, and demonstrated the following capabilities: 

1) It is much faster than traditional Monte Carlo simulation, for a comparable accuracy. 

Depending on the type of option, the computation time on a PC lies between 1 and 15 
seconds, for a 0.1% precision. 

2) It can deal with the early exercise condition of American options. To the best of our 

knowledge, this is the first method capable of pricing American path-dependent options. 

3) When the stock price 5 follows a Cox-Ross-Rubinstein binomial process, the method 

computes exactly the price of any (e.g. minimum or maximum) American lookback 
option on 5. This is also true for barrier options, such as up-and-in or down-and-out 
options. 

The next section briefly reviews the basic ingredients of the modern contingent claim valuation 
model, and discusses the alternative implementations of the numerical solutions. Section 3 
examines how path dependency can be alleviated through state augmentation, allowing for 
path-dependent claims to be priced in the standard valuation framework. It also emphasizes the 
difficulty of solving the resulting degenerate equations with standard finite difference methods. 
Sections 4 and 5 establish necessary and sufficient conditions of holonomy for degenerate dif- 
fusions. Section 6 introduces the Forward Shooting Grid method which efficiently copes with 
this degeneracy, as well as with the early exercise condition of American path-dependent secu- 
rities. Section 7 establishes the convergence of the FSG method for general multidimensional 
diffusion processes. Finally, Section 8 gives several examples of European and American 
prices for various path-dependent options, and discusses the results. 

2 Contingent claim valuation 

This section first briefly reviews the basic ingredients of the modern contingent claim valuation 
model in a continuous time framework (see e.g. Duffie (1992)), then discusses the alternative 
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implementations of the numerical solutions. 
2.1 The PDE approach 

There exists an asset-price stochastic process S t governing the evolution of the state variable 5, 
which follows the Ito's stochastic differential equation (SDE): 

dS t = m(t,S,)dt + b(t,S t )dw, S 0 > 0, (2.1) 

where w is a standard Brownian motion. 

There also exists a bond-price process B t , governed by the SDE: 

dB t = r(t, S,)B,, Bo > 0. 

This process models the continuously compounding risk-free interest rate. 

Consider a derivative security with terminal payoff g(Sr), where g is some continuous real 
function that depends on state variable 5. It is shown that in the absence of arbitrage, the price 
C(t, S t ) of the derivative security at time t solves the partial differential equation (PDE): 

BC BC 1 d 2 C 

- r(t, S) C(t, S) + —(t, S) + r(t, S) S —(t, S) + -b(t, S) 2 -^(t, 5) = 0 (2.2) 

with the boundary condition C(T, S) = g(S). 



2.2 The Feynman-Kac formula 



Under suitable technical conditions for functions r, a, and g, the Feynman-Kac (FK) formula 
gives the solution to eq. (2.2): 



C(t,S) = E, 



exp - / r(u,S u )du g(S T 



(2.3) 



where S u is the Ito process defined by: 



S u — St, u < t 
dS u = r(u, S u )S u du + b(u, S u )dw, u> t. 



Therefore, the value (price) C(t, S) of a derivative security can be interpreted as the expectation 
of its discounted payoff under the modified (aka risk-neutral) process S u , whose expected rate 
of return is the riskless interest rate of the market. 

Both the fundamental PDE (2.2) and its corresponding FK formula (2.3) extend to the case of 
a derivative security contingent to an arbitrary number of assets. 
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2.3 The Black-Scholes equation 

In the Black-Scholes model, the stock price 5 follows a log-normal process, with r and a 
constant: 

dS t = fiS t dt + aS t dw, Sq > 0 
dB t = rB t , Bo > 0. 

In this case, eq. (2.2) reduces to the Black-Scholes equation: 

f)C f)C 1 f) 2 C 

-rC(t,S) + —(t,S) + rS—(t,S) + -a 2 S 2 ^ I (t,S) = 0, (2.4) 

with boundary condition C(T, S) = g(S) = max(0, S - K). 
The change of variable 

Z = log(5/5 0 ) - at, a = r - a 2 /2, (2.5) 

simplifies eq. (2.4) into: 

dC 1 2 d 2 C „ 

- rC+ * + r^= 0 ' <2 - 6) 

with the boundary condition C(T, Z) = max(0, S 0 e z+o( - T -^ - K). 

Looking for a solution of the form C(t, Z) = e~ r ^ T ~^ E(t, Z) , eq. (2.6) reduces to the (backward) 
heat equation: 

dE 1 2d 2 E 

-m-rw = °' (2 - 7) 

with the boundary condition E(T, Z) = e r ( r " f ) max(0, S 0 e z+a ( T -') - K). 
Eq. (2.7) admits for solution the Black-Scholes option pricing formula: 

x = ^^=(io g (5/^ + (« + cT 2 )(r-o), 

where <I> is the cumulative standard normal distribution function. 

2.4 Numerical solutions 

The PDE approach and FK formula give two related ways of pricing a derivative security 
when no closed form solution exists: either solve PDE (2.2) numerically, or use Monte Carlo 
simulation (Hammersley and Hanscomb (1964)) and sample the risk-neutral process given by 
eq. (2.1) in order to compute an approximation to the expectation of eq. (2.3). The respective 
advantages and limitations of these two methods are briefly presented below. 
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Monte Carlo simulation. 

• The space complexity (memory requirement) is linear in the number of state vari- 
ables. In general, the time complexity (computation time) is quadratic in the 
number of state variables. 

• Due to the high number of samples usually required to get sufficient precision, 
execution times are significantly greater than those of finite difference methods. 

• The early exercise condition of American options cannot be dealt with. 

PDE integration. In the context of asset pricing problems, PDEs usually reduce to simple 
advection-diffusion equations, whose solutions are most simply computed through finite 
difference methods. 

• Both the time and space complexities of FD methods grow exponentially in the 
number of state variables. This makes the PDE approach attractive only when the 
number of underlying assets is small. 

• Numerical instabilities of FD methods can be a delicate issue. In particular, it is 
well-known that FD methods are ill suited to solving degenerate PDEs, that is PDEs 
for which the covariance matrix is singular. 

• In general, FD implementations run significantly faster than corresponding Monte 
Carlo simulations. 

• PDE solving methods are the only ones that easily handle the early exercise condi- 
tion of American options 1 . 

3 Path dependency 

This section first examines popular examples of path dependent problems, then describes how 
path dependency can be alleviated through state augmentation. It also emphasizes the difficulty 
of solving the resulting PDEs with standard FD methods. 

3.1 A menagerie of path-dependent options 

Historically, many similar path-dependent options have been given different names. This 
section presents the most popular path-dependent options, giving explicit formulas for their 
associated payoffs. 

Without loss of generality, we will assume all options to be issued at initial time t = 0, expiring 
at time T > 0. We will write Sj the value of asset S at expiration time T, and K the strike 
or exercise price of the option. Cj (resp. Pj) will denote the payoff of a call (resp. put) at 
expiration time, and C (resp. P) the sought value of the call (resp. put) at initial time. 

'The existence of solutions to advection diffusion PDE with early exercise boundary conditions has been studied 
by several authors McKean (1965); Merton (1973); Harrison and Kreps (1979); Bensoussan (1984); Karatzas 
(1988); Jailletetal. (1988). 
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Let M t (resp. m t ) be the maximum (resp. minimum) value of asset 5 over the time period [0, t]\ 

M t = max 0 < T < f 5 T 
m, = min 0 < T < f 5 T . 

Let also A t be the value of the arithmetic average of 5 over [0, t\: 

1 f 

A t = - S T dr (3.2) 
1 Jo 

Ordinary option. An ordinary call (resp. put) option gives its owner the right to buy (resp. 
sell) S at strike K: 

C T = max(0, S T - K) 
P T = max(0, K - S T ). 

Lookback. A lookback call (resp. put) gives its owner the right to buy (resp. sell) S at its 
lowest (resp. greatest) price over time period [0, T\: 

Cj = Sj — ynT 
Pj = Mj — Sj. 

Option on extrema. A call (resp. put) on maximum (resp. minimum) is like an ordinary 
option, but the spot is to be replaced with the historical maximum (resp. minimum) value 
of 5: 

C T = max(0,M r - K) 
P T = min(0, K — m?) . 

Capped options. A capped call (resp. put) option is like an ordinary option, as long as the 
historical maximum (resp. minimum) of S stays below (resp. above) a predefined upper 
(resp. lower) barrier price b. Should the maximum (resp. minimum) reach (resp. fall 
below) that barrier, the option is automatically exercised. 

f b-K if M r < b 
T \ max(0, St — K) otherwise 

f K-b if m T > b 

[ max(0, K — St) otherwise. 

Barrier options. Barrier options are also known as knock-out, knock-in, or trigger options. 
We describe below the most popular types. 
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• Down-and-out call. Up-and-out put. 

A down-and-out call (resp. up-and-out put) behaves like an ordinary call (resp. 
put) as long as the historical minimum (resp. maximum) of 5 stays above (resp. 
below) a predefined lower (resp. upper) barrier price b. Should the minimum (resp. 
maximum) reach or fall below (resp. reach or rise above) that barrier, the payoff 
becomes zero: 

c _ { max(0, S T - K) if m T > b 
\ 0 otherwise 




max(0,^-5 r ) if M T < b 
0 otherwise. 



• Down-and-in call. Up-and-in put. 

A down-and-in call (resp. up-and-in put) behaves like an ordinary call (resp. put) 
as long as the historical maximum (resp. minimum) of 5 stays below (resp. above) 
a predefined barrier upper price b. Should the maximum (resp. minimum) reach or 
rise above (resp. reach or fall below) that barrier, the payoff becomes zero: 




max(0,S T -K) if M T < b 

0 otherwise 

max(0, K - St) if rriT > b 

0 otherwise 



Average-rate options. Average-rate options are also known as Asian options. Their payoffs 
depend on the value of the arithmetic average of 5 over a given time period. There are 
two types of such options: fixed-strike, and floating-strike. 

• A fixed-strike average-rate option is like an ordinary option, with the time average 
At substituted for St- 

C t = max(0,A r - K) 
P T = max(0,K- A T ). 

• Symmetrically, in the payoff of a floating-strike average-rate option, the time 
average At is substituted for the strike K: 

C T = max(0, S T - A T ) 
P T = max(0,A r - S T ). 

3.2 An introductory example 

It should be noted that the PDE (2.2) for the price process can only be written down under 
the Markovian assumption that the instantaneous payoff of the derivative security be function 
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only of the current value of asset 5, i.e. C = C(S t , t). If this is not the case, the security is 
termed path-dependent. To understand the implications of such path-dependency with respect 
to the PDE solving approach, the following sections draw on an example developed in Stanton 
(1989) for the valuation of a simplified zero-strike Asian option in the Black-Scholes model. 

Consider a stock 5 with no dividends, which follows a log-normal price process. Also, assume 
a constant risk-free interest rate. An option on S is issued at time t = 0, expiring at time T > 0, 
whose terminal payoff gj is the arithmetic average Aj of 5 over period [0, T], cf. eq. (3.2). 

Applying the Feynman-Kac formula, the arbitrage-free price of this option can be computed 
as the discounted expectation of its payoff under the risk-neutral process S u : 



S t , u < t 
dS u = rS u du + aS u dw, u > t; 



and 



T 

du 



C(t) = l -e-^-')E t I S u . 
1 Uo 

= l -e-< T -<) f T Et[S u ]du 
1 Jo 



By definition of S u , E t [S" u ] = S t for u < t. For u > t, S u is a log-normal process, whose 
expectation at time t is given by: 



E t [S u ] = Ste^-'l 



Therefore, 



/ E,[S u ]du = / E t [S u ]du+ I E t [S u ]du 
Jo Jo Jt 

= ^Sudu+^e^-')- 



Finally, 



C(t) = C(t,S t ,A t ) 



-r(T-t) ft a (3.3) 

' S u du+^ l-e-^ 



o rT 

which clearly shows that the value of the security, as well as its payoff are path-dependent. 
3.3 Augmenting the state space 

State augmentation is a classical method for converting path-dependent problems into their 
equivalent path-independent counterparts. The following example is drawn from Stanton 
(1989). 
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Let us therefore incorporate A as the second state variable. A t being the average value of S over 
period [0, t],t < T, with Aq = So, the law of evolution of A t is obtained by differentiating eq. 
(3.2): 

dA t =- t (S,-A t )dt. (3.4) 

The value of the option at time t becomes C = C(S, A,t). Using the same arbitrage arguments 
as in the regular case, it is shown that in absence of arbitrage, A t solves the PDE: 

dC 0C 1 2 ,d 2 C 1 t dC n 

- rC + 7n +rS es + 7 S V ^ 7 = (3 5) 

with the boundary condition C(T, S,A) = g = A. 

Eq. (3.5) is in fact the Black-Scholes equation (2.4), into which the term l/t(S t - A t )§£ has 
been incorporated. By augmenting the state space, the value of the option only depends on 
the current values of the state variables 5 and A, hence path dependency has been alleviated. 
Furthermore, for this particular boundary condition, a closed-form formula for C(t, S t ,A t ) is 
easily found: 

C(t,S t ,A t ) = f^^A t +^(l-*-'( r -')), (3.6) 
which is exactly eq. (3.3) above. 

In the case of a fixed-strike Asian option with a more realistic terminal payoff g = max(0, A - 
K), eq. (3.5) is still valid but the new boundary condition makes it impossible to find a closed- 
form formula for the price. In general, only numerical methods can handle path-dependent 
asset pricing problems, with the usual alternative of Monte Carlo simulation vs. PDE solving. 

3.4 Degeneracy of augmented PDE 

Augmented PDE are degenerate in nature. Intuitively, the degeneracy comes from the fact 
that the augmented state variables, such as A t and S t , are correlated. In other words, the 
increments dS t and dA t of asset price and time average are not independent, which implies 
that the deterministic relationship between them, as expressed in eq. (3.5) cannot be readily 
integrated. In practice, any numerical scheme that approximates an augmented PDE, although 
it can be made stable, will not necessarily converge to a solution. This makes the standard FD 
approach impractical. Again, we shall illustrate this with the example of the zero-strike Asian 
option, eq. (3.6) providing analytic solutions for reference. 

In eq. (3.5), a discontinuity occurs at t = 0, which can be eliminated by the change of variable 
Y = tA. Together with the change indicated in eq. (2.5), they simplify eq. (3.5) into: 

dE 1 2 d 2 E BE n m 

m + 2°W + s »Y = 0 - <3 ' 7) 
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with the boundary condition E(T, Z,Y) = e Y/T. 

Because there is no second-order term in Y, eq. (3.7) is a strongly degenerate two-dimensional 
advection-diffusion PDE. invert a singular matrix. Writing E(t, Z, Y) as E" k = E(t n , Z", Y%), 
eq. (3.7) is approximated with the following explicit scheme: 

1 i 

-D,(t) = -a 2 D n (t + At)+S(t + At)D 2 (t + At), 

where 

F" _ 9F« _|_ F« fl2 

OnW - - ^ £ +°( AZ )- 

Eq. (3.7) is then backward integrated as follows: 



(7 2 Af 



" 5 i AT' 

= (l-ife^ + ^^ + ^ + ^r^-^) (3.8) 

n = N-l,...,0 

j = -n,...,n 

k = 0, . . . , k m 



with the boundary condition 

Ej,k = e- rT Y?/T 

j = -N,...,N 

k — 0, . . . , kfyi. 



We implemented the scheme (3.8), and compared results with those obtained with formula 
(3.6). Not surprisingly, the numerical procedure failed to give accurate results on a broad range 
of input parameters a and r. 

In summary, although path dependency can be alleviated by augmenting the state space, the 
applicability of FD methods in such cases is very limited. This seems to give a definite 
advantage to Monte Carlo methods, which just simulate the asset-price process. This fact 
has been noted in previous work (Talay (1991)). In Section 6, we will present an alternative 
numerical technique for solving degenerate PDEs that has several advantages over the Monte- 
Carlo method. Before, we derive a mathematical characterization of path-dependence. 
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4 Characterization of path-dependent price processes in two dimensions 



We establish necessary and sufficient conditions under which degenerate diffusion PDE can 
be reduced to lower-dimensional non-degenerate PDE on a sub-manifold of the state space. 
A degenerate PDE that can be reduced in such a way is called holonomic. It is called non- 
holonomic otherwise. It is important to determine whether a given PDE is holonomic or not, 
since this property influences both the type of numerical integration technique that can be used, 
and the memory required to conduct the integration. We introduce in this section the main 
intuition behind this characterization in the two dimensional case. In the next section, we treat 
the general case. 



4.1 The generic two-dimensional path-dependent pricing problem 



We consider a price process S t following the Ito SDE (2.1), and an arbitrary path-dependent 
variable A t : 

A t = v (rs T ] T < f ) 

We assume that A t follows an Ito SDE of the type: 

dA, = m A (t,St,A,)dt + b A (t,S t ,A t )dw 
For example, if A t is the mean value of S t up to time t, we see from the previous section that: 

m A (t,S,A) = ^(S-A), b A (t,S,A) = 0 

In order to compute the price of a path-dependent contingent claim C(f, 5 f , A t ) , we must solve 
for the following two-dimensional degenerate PDE: 



dC dC l l2 d 2 C dC l l2 d 2 C , , d 2 C 
-d7 + ms dS + 2 b *W + MA dA + 2*W + bsbA dSdA 



where: mj = rS, b$ = b. 

The numerical integration of the above equation requires a priori to quantize both variables 5 
and A. Hence, if each variable is quantized using N samples, the memory requirement of the 
numerical pricing procedure is proportional to N 2 . 

However, since the time evolutions of 5 and A depend on the same Brownian motion W, there 
is a simple relationship between the three (stochastic) differentials dt, dS, and dA: 

duj = (m s b A - b s m A )dt - b A dS + b s dA = 0 (4.2) 

The numerical integration of the PDE will require a memory proportional to N 2 only if this 
differential relationship cannot be integrated, i.e. if there is no potential function F{t, 5, A) and 
no integrating factor A(f, 5, A) such that: 

dF = Xdcj 
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Indeed, if such F and A exist, then the differential relationship (4.2) is equivalent to the 
following: 

F(t,S t ,A t )=F(0,S 0 ,A 0 ) 

Hence, we can perform (under suitable technical conditions) the following change of variables 
inPDE (4.1): 

t = ® t (t,S,A) = t 
£ = 0 s (t,S,A) = S 
F = & F (t,S,A) = F(t,S,A) 



and get: 



dC dC l l2 d 2 C dC l l2 d 2 C , , d 2 C 
d^ + ms dZ + 2 bs m + MF dF + 2 bF dF> + bsbF dWF 



where we have defined for notational convenience: 

dF dF 1 2 d 2 F dF 1 , d 2 F d 2 F 
mF= ^ + ms ^ + 2 bs ^ + mA dA + 2 bA dA^ + bsbA dSdA 

and 

dF dF 
bF = bs dS +bA dA 

But by Ito's formula: 

dF = ntpdt + bpdw 

Since by assumption dF = Xdcj = 0, we get nip = bp = 0. Hence, the transformed PDE (4.3) 
simplifies to the usual one-dimensional Black-Scholes equation: 

dC dC l l2 d 2 C n 
- rC+ lTr +mS lii + 2^m= 0 (4A> 

In order to price a contingent claim with terminal payoff g(Sr, At), one must first solve for 
A = 4>{t,S) as a function of t and 5 in equation F(t,S t ,A t ) = F(0, So,Aq), and then solve 
backwards in time the above one-dimensional PDE with the boundary condition C{T,S,A) = 
g(S,4>(T,S)). The memory required to conduct the numerical integration is now proportional 
to N instead of Af 2 . Also, the above one dimensional PDE is not degenerate, and can be 
integrated using standard FD methods. 



4.2 Fisk-Stratonovitch stochastic differentials 



It is an important matter to characterize the conditions under which the differential relationship 
(4.2) can be integrated. For that purpose, we will first introduce the notion of stochastic 
differential in the Fisk-Stratonovich sense. 

Let 0 = to < h < ■ ■ ■ < t m = t be an arbitrary partition of the interval [0, t] with a modulus 

i] = max/(f i+ i - ti). 
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Formally, the Fisk-Stratonovitch SDE dx = b o dw for two continuous semi-martingales b, w 
can be interpreted as a notation for the identity 



m— 1 



V/, 



x(t) = *(0) + lim V %)+ 2 %+l) (^ + i) " 

v HI * — ' Z 



where the limit is taken in probability. 

Intuitively, whereas the Ito differential dx = bdw is a. forward differential which can be loosely 
interpreted as: 

x(t + dt) = x(t) + b(t) (w(t + dt) - w(t)) 

the Fisk-Stratonovitch differential dx = b o dw is a symmetric differential, which can be loosely 
interpreted as: 

bit) + bit + dt) . . 
x(t + dt) =x(t) + ^ ^(w(f + <*) - w(t)) 

It can be shown (see e.g. Karatzas and Shreve (1988)) that, unlike Ito differentials which 
follow Ito's chain rule, Fisk-Stratonovitch differentials follows the chain rule of the classical 
differential calculus. In particular, for any smooth real-valued function / (t, x\, . . .,Xd) and 
continuous vectors of semi-martingales X = (jci , . . . , xj) , we have: 

df(h X) = f t (t, + £ J£ (t, X) o dx t 

i=\ 1 

Finally, it can be shown that there is the following relationship between Ito and Fisk- 
Stratonovitch SDE. If X = (xi, . . . , x d ) follows the Ito SDE: 

k 

Vi G [l,d], dxi = nti(t,X)dt+ y^ j b i j(t,X)d\Vj 

Then 

d k d£> \ * 
m i ~ 2 zC zC &f ] * + zC ^ ° ^ 

1=1 j=l J y=l 

4.3 The stochastic Frobenius integrability condition in two dimensions 

In this subsection, we introduce the integrability conditions for equation (4.2). We consider 
for simplicity the case d = 2, k = 1. The system of Ito SDE: 

dx\ = m\ it, x\,X2)dt + b\ it, x\,x-i)dw 
dx2 = m-iit, x\,X2)dt + b-iit, x\,xi)dw 

admits the following Fisk-Stratonovitch equivalent: 

dx\ = m\ it, x\,x-i)dt + b\it, x\,xi) o dw 
dx2 = m-iit, x\,x-i)dt + biit, x\,xi) o 
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with 



m 2 = m 2 -\(b 1 ^ + b 2 ^ 



2 dxj ~ " z dx 2/ 

By eliminating dw, we get the following differential relationship: 

dcj = cjQdt + wio dx\ + lo 2 o c?X2 = (mi^2 - m 2 b\)dt - b 2 o dx\ + b\ o dx 2 = 0 

By the Fisk-Stratonovitch (i.e. classical) chain rule, a sufficient condition for the existence of 
a potential function F(t, xi,x 2 ) and an integrating factor X(t, xi,x 2 ) such that dF = Xdcj is: 

dF dF dF 

— = \u 0 , — = Xcji, — = Xlj 2 (4.5) 
at ox\ ox 2 

From the above relations and the symmetry of the second order mixed partial derivatives of F, 
we get after elimination of A: 

Uo ( —1 - dtjJl \ ^1 - +u 2 { — - dtjJ ° \ - 0 (4 6) 

\ dx\ dx 2 ) \ dt dx 2 ) \ dt dx\ ) 

It can be shown (see section 5) that the above relation is a necessary and sufficient condition for 
the existence of F. We can express the above condition in terms of m; and bj instead of a;,-. Let 
us denote Bo = (1> m 2 ) T and Bi = (0, bi,b 2 ) T . After elementary algebraic manipulations, 
the integrability condition (4.6) is shown to be equivalent to: 

det (£(,,#!, [So, #i]) =0 

where [Bo, #1] is the Lie Bracket (see also section 5) of the vector fields Bq and B\, i.e. the 
vector field: 

[BqiBi] = dBo.Bi - dBi.Bo 
and where dBj denote the Jacobian matrix of the vector field Bj. 



4.4 Characterization of path dependent prices processes 



We consider the path-dependent price process of subsection (4.1). We take x\ = 5 and x 2 = A. 
Define 

1, db s 

m s = m s - -b s — 

1 (u db * , u db A 
mA = mA -2\ bs ^S +bA M) 

and let Bq = (1, %, fhp) T , and B\ = (0, bs, bA) T ■ Then: 



The degenerate two-dimensional Black-Scholes PDF (4.1) in the variables (S,A) 
can be reduced to the one-dimensional non-degenerate PDF (4.4) if and only if 
the following condition is satisfied: 

det (B 0 , S l5 [B 0 , Si]) =0 
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We can apply the above result to the case of a path-dependent variable A which is instantaneously 
riskless. i.e. b A = 0. 

We have: Bo = (1, rhs, «u) T , B\ = (0, bs, 0) T , hence 



Therefore, the problem can be reduced to a one-dimensional problem iff m A does not depend 
on 5. This is obviously not the case for the average-rate option pricing problem: 

dmA = 1 

ds t T 

We can state: 

The average-rate option pricing PDE (3.5) cannot be reduced to a one-dimensional 
non-degenerate PDE. 

5 Characterization of path-dependent price processes: general case 

In this section, we turn to a generalization of the above result. 

5.1 The general diffusion problem 

We consider an advection-diffusion equation of the type: 



where X = (jci , . . . , x d ) T G R d , M = (mi , . . . , m d ) T , with the following boundary condition: 

f(T,x h ...,x d ) = g(xi, ...,x d ) 

We assume that the real valued functions r(t,x\, . . . ,x d ) > 0, m,(f,xi, . . . , x d ) , i £ [1, d], and 
jij(t, x\ , . . . , x d ) , (ij) G [1 , d] 2 on [0, oc[xR d are bounded and sufficiently smooth. We assume 
furthermore that for any t, x\ , . . . , x d the matrix: 



det (B 0 , B h [B 0 , Si]) = -b] 



2 dm A 



s dS 




(5.1) 



Y(t,x\, ...,x d ) = (jij(t,xi, . . •,^))( ; - l /)e[i,^2 



is symmetric and non-negative. We denote by B its Cholesky decomposition: 



r = bb t 



B is a d X k matrix of rank k, with k < d. 
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By the Feynman-Kac formula, the solution of the above PDE can be expressed: 
f(t,x h ...,x d )=E t (exp(- J r(T,X(r))dT^g(X(T)) 

where X is the solution of the SDE: 

dX = Mdt + BdW, X(t) = (x u ...,x d ) T (5.2) 

W = (w\ , . . . , Wk) is a ^-dimensional standard Brownian motion. 
We can define the vector M: 

1 ^ ^ db 
Vi G [1, d], m = tm - - V V by— ^ 

Then, from subsection (4.2): 

tff = + BodW (5.3) 

We denote by S w (Xq) the support of the above SDE for X(0) = Xq, i.e. the smallest closed set 
of paths X(t) (elements of the Wiener space Karatzas and Shreve (1988); Ikeda and Watanabe 
(1981)) such that ¥xob(S w (Xq) ) = 1. The set of accessible states at time t of SDE (5.3) for 
X(0) =Xo,i.e. the set of values X(t) for elements of S w (Xo), is denoted by A w (t, Xo). In other 
words, A w (t, Xq) is the time t projection of S w (Xq). Clearly, 

Prob(X(f) G A w (t,Xo)) = 1 

The set of accessible states from Xo for system (5.3) is defined as: 

A w (X 0 ) = U t >oA w (t,X 0 ) 

Clearly, 

Prob(V t,X(t) eA w (X 0 )) = 1 

5.2 The certainty equivalence theorem of Stroock and Varadhan 

The following theorem was first established by Stroock and Varadhan (1972). We use here 
the formulation of Ikeda and Watanabe, where the theorem is established as a consequence of 
general convergence results for approximations of diffusion processes (see Ikeda and Watanabe 
(1981)). 

We consider the deterministic control system: 

dX - dU 

— = M + B— (54) 
where U (t) = (u\(t) , . . . , Uk(tj) is a piecewise smooth function of time. 
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We define the support S u (Xq) of system (5.4) as the set of all possible paths X(t) starting at 
Xq for all possible control paths U. Similarly, we define the set A u (t, Xq) C R d of accessible 
states at time t from Xq as the set of states X* £ R d such that there exist a control U (t) 
verifying X* = X(t). The set A u (Xq) c R d of accessible states from Xq is similarly defined 
as A u (Xq) = U f >oA f/ (f,Xo). Finally, we say that a state X* is weakly accessible from Xo 
if there exist a sequence Xo,Xi, . . . , X n = X* such that either is accessible from Xi, or 
Xi is accessible from The set of states weakly accessible fromXo is denoted WA(Xq). 

By definition, any accessible state is weakly accessible, i.e. A(Xq) c WA(Xq). For a non- 
symmetric system, i.e. a system such that the controls U cannot always be inverted, then there 
may be weakly accessible states that are not accessible. 

We have the following fundamental result (Stroock and Varadhan (1972)): 
(Stroock and Varadhan, 1972) 

For any Xq, the support of the Fisk-Stratonovitch system (5.3) is the closure of 
support of the deterministic system (5.4): 



We will only use the following immediate corollary: 

For any Xq, and any time t, the set of accessible states at time t from Xq of the 
Fisk-Stratonovitch system (5.3) is the closure of the set of accessible states at time 
tfrom Xq of the deterministic system (5.4): 



5.3 Frobenius integrability condition and Chow's theory 

In the previous section, we used the certainty equivalence theorem in order to establish a link 
between stochastic and deterministic accessibility for differential equations. In this section, 
we apply standard result on the accessibility of deterministic control systems to characterize 
the integrability of degenerate SDE. 

We consider again the deterministic system (5.4). For the purpose of the following discussion, 
we will include the time variable t in the state, i.e. we consider the state of the system at time t 
to be Y = (yo(t),yi(t), ...,yd(t)) = (t,xi(t), . . .,x d (t)) e R d+l . If by denote the elements of 
the matrix B, we define the extended (d + 1) X (d + 1) matrix B: 



(Xq) = SU(X 0 ) 



W>0, VX 0 , A w (t,X 0 ) =AU(t,X 0 ) 



( 1 ° 

rh\ b\\ 



0 \ 



B = 



\ m d b d \ 



bdd ) 
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Then, system (5.4) is equivalent to the following system in R d+l : 

dY ~dU 

n,= B m <5 - 5) 

where U(t) = (t,U(t)) = (t, u\ (t),..., Uk{t)) = (&o(f), ui(t), . . . , is the extended 

control. 

Let A u (Yq) denote the set of extended states accessible from Yq = (0, Xq). Clearly: 

A°{Y Q ) =U ( > 0 {/} xA u (t,Xo) 

and 

WA D (Y 0 ) = U f >o{0 x WA u (t,X 0 ) 

By definition, any accessible state is weakly accessible. Reciprocally, if the drift is zero, i.e. 
M = (mi, . . . , fhd) = 0, then weak accessibility is equivalent to accessibility. Indeed, since 
the system (5.5) is symmetric in the control variable U (i.e. — U is an admissible control iff U 
is an admissible control), the extended control U becomes symmetric for a zero drift. 

If M = 0, we get: 

W, WA u (t,Y 0 ) =A u (t,X 0 ) 

We now recall the definition of the Lie bracket of two vector fields. Let (V 1 , V 2 ) be any pair 
of vector fields in R d+l . 

Given any point Yo = (t,Xo) G R d+l , let us consider a path starting at Yo and obtained by 
concatenating the four following paths: 

- the first path follows the flow 2 of V 1 during St; 

- the second path follows the flow of V 2 during St; 

- the third path follows the flow of — V 1 during St; 

- the fourth path follows the flow of -V 2 during St. 

Let Y\ = (t + St, X\) be the configuration reached at the end of these four paths. A straightfor- 
ward Taylor expansion shows that: 

Km Xizll = dV 2 • v l - dV l ■ V 2 , 

St^Q St 2 

where dV 2 ■ V 1 and dV 1 ■ V 2 denote the products of the (d + 1) X (d + 1) Jacobian matrices: 



2 The integral curve of a vector field V is a curve whose tangent at every point Y is V( Y). We say that the curve 
follows the flow of V. 
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dV 2 



H 

dyo 



9v| 
dyo 



H 
dya 



9v| 
dyd 



dV 1 



dyo 



\ dy 0 



dy d X 



wy d 
dyd 



and the (d + 1) -vectors: 

V' = {vlv\ ... vtf; V 2 =(v 2 v 2 ... v 2 ) T . 



The expression dV 2 ■ V 1 — dV 1 ■ V 2 determines a new vector field which is commonly denoted 
by [V 1 , V 2 ] and called the Lie bracket of V 1 and V 2 . 

Let Bq, . . . , Bd be the columns of the extended matrix B. By definition, the Control Lie Algebra 
associated with system (5.5), denoted by CLA(B), is the smallest Lie algebra which contains 
Bq, . . .,Bd- Stated otherwise, CLA(B) is the subspace of vector fields generated by all the 
linear combinations of vector fields Bo, . . . , Bj and all their Lie brackets recursively computed. 

For every Yq G R d+l , let CLA(B)(Yq) denote the subspace of vectors spanned by the vector 
fields of CLA(B) at Yq. A connected sub-manifold M. of R d+l is an integral sub-manifold 
of CLA(B) if at each Y £ R d+l the tangent space to M is contained in CLA(B)(Y). M is a 
maximal integral sub-manifold of CLA(B) if it is not properly included in any other integral 
manifold. 

The Frobenius integrability theorem can be stated as follows: 

Frobenius Integrability Theorem If the dimension ofCLA(B) (Y) has a constant 
value rfor every Y £ R d+l , there exists a partition ofR d+l into maximal integral 
sub-manifolds of CLA(B) all of dimension r. 

The maximal integral sub-manifold of dimension r passing through Y is denoted M r {Y). 

The following results derive from the work of Chow (1939). They were elucidated in Hermann 
(1963); Haynes and Hermes (1970); Lobry (1970); Sussmann and Jurdjevic (1972); Hermann 
and Krener (1977). We follow the presentation of Hermann and Krener (1977). 

Chow's Theorem If the dimension of CLA(B) (Y) has a constant value r for every 
Y e R d+l , then 

WeR d+ \ Mr{Y) = WA U (Y) 
Furthermore, the interior of A u (Y) as a subset of M. r (Y) = WA U (Y) isnotempty. 
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5.4 A characterization of holonomy for degenerate diffusions 

By inspection of the extended matrix B, and since the original matrix B is by definition regular, 
the minimal dimension of the Control Lie Algebra is r m; „ = k + 1 . 

Let us assume that the Control Lie Algebra has a constant dimension k+l<r<d+l. 
We will show that there exist exactly m = d + 1 - r constraints satisfied almost surely by the 
solution of the Fisk-Stratonovitch system (5.3). 

We will first show m > d + 1 - r. Indeed, Chow's Theorem states that the set of weakly 
accessible states for (5.5) is a sub-manifold of dimension r. This means that there exist a set of 
d+ 1 - r independent constraints F = [Fi(t,xi, . . .,Xd), . . . , Fj + i- r (t,xi, . . .,Xd)) such that 

VteR, vx 0 eR d , WA u (t,x 0 ) = {xeR d , F(t,x) = f(o,x 0 )} 

This result, together with the Certainty Equivalence Theorem, implies that the set of accessible 
states fromXo of the Fisk-Stratonovitch system (5.3) satisfies the constraints F(t, X) = F(0, X$) 
almost surely. Indeed, the set of accessible states is a subset of the set of weakly accessible 
states. Hence m > d + 1 — r. 

Reciprocally, let us assume that the solutionX(f) of the Fisk-Stratonovitch system (5.3) satisfies 
almost surely m independent constraints F = (Fi, . . . ,F m ). By the certainty equivalence 
theorem, we have F(A u (Yq)) = F(Yq). But, by Chow's theorem, the interior of A v (Yq) in 
M. r is not empty. Hence, it cannot satisfy more than d + 1 — r constraints, since M. r is of 
dimension r. This implies m < d + 1 — r. 

We can state: 

Condition for the existence of d + 1 - r constraints 

If the Control Lie Algebra CLA(B) has constant dimension r, then there exist 
exactly d + 1 - r independent constraints F = (Fi, . . . , F^ + i_ r ) such that the 
solution of the Fisk Stratonovitch system (5.3) satisfies the constraints almost 
surely. 

Prob (Vf,F(f,X(f)) = F(0,X 0 )) = 1 

A degenerate Ito SDE of the form (5.2) and a degenerate advection-diffusion PDE on the 
form (5.1) are called holonomic if they can be reduced to lower dimensional non-degenerate 
counterparts on an appropriate sub-manifold of R d . They are called non-holonomic otherwise. 
We have the following characterization of holonomy: 

Characterization of holonomy for degenerate diffusions 

Consider a degenerate d-dimensional ltd SDE of rank k < d of the form (5.2) and 
its corresponding degenerate advection-diffusion PDE on the form (5.1). Assume 
that the control Lie algebra CLA(B) has a constant dimension r. 

Then, the ltd SDE (5.2) and the corresponding PDE (5.1) are holonomic if and 
only ifr = k+\ 
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When dealing with a degenerate PDE of the form (5.1), we must first compute the dimension 
of the control Lie algebra. If the dimension is minimal (k + 1), then we must change variables, 
replacing the last d — k variables by the d—k integrable constraints F\ , . . . , F^-k- Then, the PDE 
is transformed into a non-degenerate PDE on the k dimensional integral sub-manifold satisfying 
the constraints. Hence, standard FD methods can be used for the numerical integration. 

On the other hand, if the dimension of the control Lie algebra is not minimal (r > k+ 1 ), then the 
PDE is non-holonomic, i.e. intrinsically degenerate. Then, FD methods are not appropriate. 
In the next section, we describe an alternative numerical integration technique for solving 
non-holonomic diffusion PDE. We first describe the method for typical two-dimensional path- 
dependent pricing problems such as the average-rate option pricing problem. Then, we present 
the method in full generality. 

6 The Forward Shooting Grid (FSG) method 

This section introduces a new numerical method for the path-dependent asset pricing problem, 
which we call the Forward Shooting Grid method (FSG). FSG efficiently copes with the 
degeneracy of augmented PDEs. More generally, the FSG method is adequate for solving any 
non-holonomic advection-diffusion equation. 

The FSG method was first introduced (under a different name) by Barraquand and Latombe 
(1993) for solving non-holonomic Hamilton-Jacobi-Bellman equations in a deterministic set- 
ting. However, to the best of our knowledge, the FSG method has never been used before for 
solving non-holonomic stochastic optimal control problems such as the advection-diffusion 
problem described here. 

6.1 Principle 

The FSG method consists in taking advantage of the forward SDE equations that govern the 
correlated evolutions of the augmented state variables with respect to the underlying asset 
variables. Combined with an a priori quantization of the augmented state, those equations 
allow to construct the discrete state graph of the correlated variables, which is used in turn to 
integrate backwards in time the degenerated equation in the state variables. 

Consider a path-dependent contingent claim on an asset 5, with terminal payoff Cj = g (St , Aj) , 
where A is a path-dependent variable, e.g. the historical minimum or average of S. Also, assume 
5 follows an Ito SDE, such as in eq. (2.1). 

By augmenting the state space with A, the price C of the claim depends on both 5 and A. Let 
us define an a priori quantization as follows. Given a time step At, we fix two values AT, AZ, 
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and two invertible quantization functions S, A, such that: 

S t = S(nAtJAZ) = SJ 

A t = A{nAl, kAY) = A\ 

(6.1) 

n = 0,...,N=T/At 

j = jo(n),...j m (n) 
k = ko(n), . . .k m (n). 

Second, let us assume a function tp quantifying the correlated evolutions of A with respect to 
5, i.e. for an arbitrary time step At: 

). (6.2) 

Eq. (6.2) relates the variations of 5 and A upon a transition from state (S t ,A t ) to state 
(S t+ At,A t+At ). Under the quantization (6.1), it has the following discrete equivalent: 

A"+ l = (6.3) 

Knew ' Y ™ Jnew J 

where j new (resp. k new ) denote the set of all j (resp. k) values used to approximate state 
(S t+ At,A t+ At) of the transition. By assumption, the law of evolution of 5 is known, hence 
jne W values can be obtained from the approximation of the corresponding SDE (e.g. through a 
Cox, Ross, Rubinstein (CRR) binomial approximation). In order to get the corresponding k new , 
we invert the quantization function A at time (n + I) At, and take the resulting nearest integer 
value: 



knew = nearest 



«+i 



Ay 



(6.4) 



= 4>{k,jnew)- 

In other words, we obtain k new by shooting the best approximating A bucket forward in time 
through function <f>, hence the name of the method. 

The last step consists in finding the law of evolution for C" k = C(nAt, S", AD . By the Feynman- 
Kac formula, the price C(0, 5, A) can be computed as the discounted expectation of its payoff 
future under the risk neutral process (S Ul A u ). Choosing a CRR binomial approximation, 
we have j new = {j + l,j — 1}. We let k new = {k+,k-}, where k- = <f)(k,j — 1), and 
k+ = 4>(k,j+ 1). This gives: 

q k = uC»+l k+ + (l-u)Cpl k _, (6.5) 
where u = 1/2 + oiAtjla 1 is the associated risk-neutral probability. 

The pricing algorithm proceeds in two steps by first building the discrete state lattice induced 
by the binomial approximation, using the forward integration of eq. (6.3), then computing the 
price C, using the backward integration of eq. (6.5). It is straightforward to check for the early 
exercise condition at each time step n of this scheme. 
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6.2 Application 

This section details the numerical valuation of two typical path-dependent instruments: a 
lookback call, and an average-rate call options on stock S. Thereafter, pricing algorithms for 
other path-dependent instruments based on the historical minimum (maximum) or time average 
of S are easily derived. This includes all types of barrier options on stock, futures or foreign 
exchange. 

6.2. 1 Lookback option 

We consider a lookback call option on a stock 5, assumed to follow a log-normal process. The 
call has terminal payoff Cj = Sj - mj, where m t is the historical minimum of S, cf. eq. (3.1). 

Given a time step At, we fix two values AZ, AY: 

AZ = aVAt 



and quantize 5 and m as follows: 



AY = AZ; 



S] = So** 2 
m n k = S 0 e kAY . 



(6-6) 



The choice AZ = AY is justified by the fact that since m is minimum over S, it necessarily 
yields one of those values. Therefore, m will be optimally quantized by choosing AY = AZ. 

The correlated evolutions of S and m arise from the definition of the minimum: 

m t+At = min(m f , 5 f+Af ) • (6-7) 

Under a binomial approximation of S and m, we associate to the upward (resp. downward) 
transition SJ -> 5"+/ (resp. SJ -> 5"+/) in 5, the transition m\ -> m n k + [ (resp. m\ -> m^ 1 ) in 
m. The discrete equivalent of eq. (6.7) is then: 

m"+ l = minK,^ 1 ), 

from which the following values k± of k+ and k— are found, using definition (6.6) of the 
quantization: 

k± = min(it,7± 1). (6.8) 

Taking the risk-neutral probability u from a CRR approximation of the log-normal process 
governing 5, the lookback call price equation writes: 

q k = uq+l k+ + (i-u)cj+l k _ 

n = N-1,...,0 < 6 " 9 ) 
j,k = n; 
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with boundary condition 

Cf k = Sf-mf = Soe^-Sof" 
j,k = -N,...,N. 

The range of variations of k± is trivially bounded, i.e. from inspection of eq. (6.8), \k±\ < n. 
The scheme (6.9) is therefore feasible. Using eq. (6.8) to compute k± values, eq. (6.9) is 
eventually backward integrated in order to get the price Cq 0 of the lookback call. 

The pricing algorithm of a lookback put is similar. The terminal payoff changes to Pj = 
Mj - St, where M t is the historical maximum of 5, cf. eq. (3.1). In this case, the quantization 
(6.6), and the recursion equation (6.8) still hold, with Mf and max substituted for m£ and min 
respectively. The price equation (6.9) becomes then: 

p- k = + (i - -WV- 

n = N-1,...,0 ( 6 - 1Q ) 

j,k = -n, . . . , n; 

with boundary condition 

Pf k = M N k -Sf = S 0 e kAY - So^ 

j,k = -N,...,N. 

6.2.2 Average-rate (Asian) option 

We consider an average-rate (Asian) call option on a stock 5 assumed to follow a log-normal 
process. We shall develop the example of a floating-strike average-rate call, with terminal 
payoff Cj = max(0, Sj - At), where A t is the time average S, given by eq. (3.2). 

Given a time step At, we fix two values AZ, A7: 

AZ = <jVa1 



AY = pAZ, p < 1; 



and quantize 5 and A as follows: 



A" k = S 0 e k ' 



„kAY 



(6.11) 



Notice that since A is an average, it does not necessarily yields one of the S values. In order 
to preserve accuracy, the quantization step AY has to be smaller than AZ, hence p < 1 . We 
postpone the determination of p until the price equation is established. 

The relation giving the correlated evolutions of A with respect to 5 arise from eq. (3.4): 

_ (t + At) A t + At S t+At 
A t+At - J ^ Kt (6.12) 
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Under a binomial approximation of S and A, we associate to the upward (resp. downward) 
transition S] -> SJ+I (resp. S] -> 5J+/) in 5 the transition A" k -> A n k + l (resp. A\ -> A n +_ l ) in A. 
The discrete equivalent of eq. (6.12) is then, with A{j = 



(n+l)A" k + Sj+l 
n + 2 

{n+l)A n + S n+l 

n + 2 



from which the following values of k+ and k— are found: 

log 

k± = nearest 



n + 2 



P AZ 



(6.13) 



Taking the risk-neutral probability u from a CRR approximation of the log-normal process 
governing 5, the average-rate price equation writes: 



q k = ucj+l k+ + (i-u)cj+l k _ 

n = N-l,...,0 

- —n, . . . , n 

'■ ~ km (w) , . . . , k m (n) ; 



(6.14) 



with boundary condition 

rN - 



5f-Af 



S 0 e>^ - S 0 e 



kAZ 



J 
k 



-N,...,N 
-k m (N),...,k m (N). 



In order for the previous scheme to be feasible, the following three points need to be addressed. 



• Bound the range of variations ofk( i.e. find k m ). 

For any time step n, the maximum value max^ A" k for average A cannot be greater than the 
maximum value max, SJ for asset S, hence A n km < S n n . Eq. (6.1 1) implies then k m < nj p. 
Thus, we let: 

* - - 

P 

• Bound the range of variations ofk+. 

Starting from eq. (6.13), and after some algebraic manipulations, one gets \k+\ < 
(n/p) + l/p, hence: 



\k+\ < k m + 



1 
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• Find p value. 

Because of the very simple structure of eq. (6.14), the propagation of error due to the 
quantization (6.11) can be shown (see next section) to be less than Soe km pAZ(e pAZ — 1). 
With AZ = a^/At, this gives a maximum relative error e ~ pAZe Na ^. Typically, 
N = 100, At = 0.00274 (1 day), a < 1, hence e ~ 0.05 p. Thus, the value p = 0.02 
ensures a 0.1 % precision. In practice, cf. sec. 8.2, the value p = 0.1 is sufficient. 



7 Convergence of the FSG method 

7.1 Lipschitz conditions 

In the previous section, we introduced the FSG method for typical path-dependent asset pricing 
problems. However, the principle underlying the FSG method is very general, and can in fact 
be applied to arbitrary diffusion equations. We analyse below the convergence of the FSG 
method for diffusions equations of the type (5.1). 

In order to prove the convergence of the FSG method, we will assume that the boundary 
condition g in problem (5.1) is Xg-Lipschitz, i.e.: 

VXeR d ,YeR d , \g(X) - g(Y) \ < K g \\X - YWoo 

with 

Halloo = max \xi\ 

ie[l,d] 

This assumption is very reasonable for option pricing applications, since most payoff functions 
are clearly 1 -Lipschitz. 

For the quantization of the PDE (5.1), we select d quantization functions whose inverse map 
the original space variables x\ , . . . , x^ into the transformed spaces variables yi, . . .,y d . 

X{t, Y)=(Xi{t,y<),...,X d {t,y d )) 

We assume that Xj is a strictly monotonic function of its second variable yu and we define its 
inverse by the following relation: 

Xr\t 1 X i (t 1 y i ))=y i 
Finally, we assume that Xj is ^-Lipschitz in its second variable y,-, i.e.: 
V/, yj, yl \Xi(t, y\) - Xi(t, yf) \ < K^yj - yj\ 

7.2 Quantization of time 

It is a well-known consequence of the central limit theorem that the standard Brownian motion 
is the limit when At — > 0 of the binomial distribution of step \J~At (see e.g. Cox and Rubinstein 
(1985); Duffie (1988)). We define the ^-dimensional binomial process W At as follows: 

Ve = . . . , e*) G {-1, l}\ Prob(W Af (; + At) = W At (t) + eVAt) = 1 
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In light of the previous discussion, we get under suitable technical conditions (see e.g. Duffie 
(1988)): 

f(t,x h ...,Xd)= JimJ(t,xi ...,Xd)=E, ^exp ^- J r(r,l(r))^ g(X(T))^j 

where X is defined by the forward stochastic difference equation: 

X(t + At) = X(t) + M(t, X(t))At + V(t, X(t))(W(t + At) - W(t)) 

Remark: The definition of X above corresponds to choosing an Euler scheme for the quanti- 
zation of the SDE. Other schemes of higher order of convergence in At could be used, such as 
Talay-Milshtein schemes (Milshtein (1974); Talay (1984)). However, we did not explore this 
alternative. 

From the law of iterated expectations, we see that/ satisfies the following recursive backward 
equation: 

~ f{tlX)= \+WY)¥ S f(t + At,X + M(t,X)At+eV(t,X)VAt) (7.1) 



where X) is defined by: 

1 



l+R(t,X) 

7.3 Quantization of space 



exp (-r(t,X)At) 



We select d quantization steps A7 = (Ay\ , . . . , Ayd) for the transformed space variables 
yi,...,yd- We denote by J the <i-uple of integers J = (ji, . . . Jd), and by JAY the <i-uple 
JAY = (j l Ay l ,...J d Ayd). 

For any given couple (n, J), and any e e {-1, l} k , we define J new (n, J, e) = (/f w , . . . ,j n / w ) 
by: 



Vi G [l,d\, j" ew = nearest 
with 



Av,- 



X* = X(nAt, JAY) + M(nAt, X(nAt, JAY)) At + eV(nAf, X(nAf, J AY)) V At 

In other words, we obtain J new by shooting forward in time the best approximating bucket in 
the space variables X. 

Then, we approximate equation (7.1) above by: 
~f a PP rox {nAt, JAY) ' 



1 +R(nAt,X(nAt,JAY)) 
T*Eee{-i,i}«f aPPr ° X (( n + ^t,X((n+ l)At,J new (n,J,e)AY)) (7.2) 
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7.4 Convergence 

Combining (7.2) with (7.1) and with the Lipschitz properties of functions g and X, we see that 
at time (N — I ) At, we have for all possible X: 

]f a PP rox ((N - l)Af,X) -f((N - l)Af,X)| < A» max K t max Ay,- 

" ' ie[l,d] ie[l,d] 

Then, a straightforward backward induction on n shows that the same inequality is true for all 
times prior to N. Finally: 

Mn < N, VX, \f approx (nAt, X) -f(nAt, X) | < K g max K t max Ay,- 

ie[l,d] ie[l,d] 

Hence, we can guarantee a precision on the result of Af by choosing the space quantization 
steps such that: 

Vje[M, Ayj<- ^ 

K g max ie[14] Ki 

If the payoff is that of an option, we have typically K g = 1. If we furthermore choose a 
logarithmic quantization, i.e. if y,-) = e yi then Kj is simply the maximum value taken by 
the variable x, at any point in time. Therefore, it is practically easy to choose the quantization 
steps Ay,- for reaching a prespecified desired accuracy. 

Furthermore, the above inequality shows that the FSG method is unconditionally convergent. 
Indeed, the approximate value converges towards the theoretical value whenever At and Ay,- 
converge towards zero. This is true regardless of any quantitative relationships between At and 
Ay,-. 

If k = d, i.e. if the covariance matrix F is regular, then finite difference methods are also 
convergent, and are faster than the FSG method. However, as soon as k < d, i.e. when 
the covariance matrix is degenerate, finite difference methods introduce a spurious numerical 
diffusion, whereas the efficiency of the FSG method is unaffected. 

The FSG method therefore is an attractive solution for all degenerate numerical valuation 
problems. Such problems arise in many other cases than the path-dependent case. 



8 Results 

This section reports European and American prices found with the FSG method for different 
path-dependent options. Since the FSG method reduces to the classical CRR method for 
ordinary options, prices of at-the-money ordinary options are also given for reference. 

The stock price So =100 and the interest rate r = 10% are held constant, while the values of 
the volatility a, the maturity T, and the exercise price K vary. Maturities 3, 6, and 12 must be 
read in months, corresponding to 91, 182, and 364 days respectively. Both call (C) and put (P) 
prices are reported. When possible, we also indicate the corresponding analytic price C a or P a . 
In all cases, we have chosen the number of time steps so as to reach a 0.1% accuracy. 
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Ordinary option 



At-the-money ordinary call option 
S 0 = 100, r = 10% 





European 


(J 


T 


analytic FSG = CRR 

C a C 






10 


3 


3.438 3.435 


6 


5.837 5.835 


12 


10.284 10.283 






20 


3 


5.286 5.280 


6 


8.262 8.258 


12 


13.244 13.241 






40 


3 


9.148 9.138 


6 


13.559 13.551 


12 


20.285 20.280 



At-the-money ordinary put option 
So = 100, r = 10% 





European 


American 


(J 


T 


analytic FSG = CRR 

Pa P 






10 


3 


0.976 0.973 


1.223 


6 


0.973 0.971 


1.450 


12 


0.792 0.791 


1.632 






20 


3 


2.823 2.818 


3.064 


6 


3.398 3.394 


3.913 


12 


3.752 3.749 


4.811 






40 


3 


6.686 6.675 


6.909 


6 


8.694 8.687 


9.205 


12 


10.794 10.788 


11.944 



Table 1: Reference table for an at-the-money ordinary option. Time step was set to 0.5 day for 
all maturities. Precision for CRR prices is 0.1%. 

Prices for at-the-money ordinary options are listed in table 1 . The time step was chosen constant 
(0.5 day) for all maturities. For moderate values of the volatility (< 50%), about N = 100 
time steps are enough to get 0.1% precision. The time complexity is 0(N 2 ), and the memory 
requirement O(N). The computation time for one call/put is a few milliseconds on a DEC 
alpha PC, and the required memory 20 Kilobytes. 

8.1 Lookback option 

Results for a lookback option are presented in table 2. The time step was chosen constant 
(1 day) for all maturities. In practice, N = 100 time steps are sufficient. The time complexity 
is 0(N 3 ), the memory requirement is 0(N 2 ). The computation time for one call/put is about 1 
second on a DEC alpha PC, and the required memory 400 Kilobytes. 

Lookback payoffs depend on extreme values of the underlying asset S. The continuous time 
framework, in which closed-form formulas are derived, captures the variations of those extrema 
over infinitely small time periods. In practice however, the extrema values in a lookback contract 
are to be computed on a daily basis. Therefore, lookback prices computed from a closed-form 
formula incorporate spurious intra-day variations, and thus overestimate the prices. Compared 
to the FSG method with a 1 day time step, this overestimation is about 5% 3 . 

3 By contrast, both the CRR approximation to the ordinary option, and the FSG approximation to the average-rate 
with 100 time steps give European prices close to within 0. 1% to their respective analytic (continuous) prices. 
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Lookback call option 
So = 100, r = 10% 





European 


American 


a 


T 


analytic FSG 

C a C 






10 


3 


5.266 5.028 


5.028 


6 


8.253 8.019 


8.019 


12 


13.298 13.075 


13.075 






20 


3 


8.939 8.482 


8.482 


6 


13.174 12.733 


12.733 


12 


19.614 19.203 


19.203 






40 


3 


16.052 15.209 


15.209 


6 


22.658 21.872 


21.872 


12 


31.804 31.106 


31.106 



Lookback put option 
S 0 = 100, r = 10% 





European 


American 


a 


T 


analytic FSG 

Pa P 






10 


3 


2.927 2.668 


3.087 


6 


3.632 3.368 


4.209 


12 


4.281 4.012 


5.630 






20 


3 


6.969 6.434 


6.845 


6 


9.283 8.729 


9.606 


12 


12.021 11.448 


13.307 






40 


3 


15.559 14.406 


14.811 


6 


21.685 20.456 


21.362 


12 


29.906 28.583 


30.661 



Table 2: Results of the FSG method for a lookback option. Time step was set to 1 day for all 
maturities. Precision for FSG prices is 0.1%. Lookback analytic prices are overestimated by 
about 5%. 



On the other hand, we checked that the above discrepancy would vanish as the time step 
tends towards 0. Results are presented in table 3 below for a 3 months lookback call option. 
Interestingly enough, while the decrease of the precision AC/ C is linear in the time step for an 
ordinary option as well as an average-rate option, it is only square root for the lookback option. 



Finally, in order to estimate the accuracy of the FSG method, we have checked by an additional 
convergence test that 1 day time step yields a least 0.1% accuracy in the computed price. 

Lookback options are more expensive than ordinary at-the-money options, roughly twice as 
much. As pointed out in Conze and Viswanathan (1991), the values of an American and a 
European lookback call are equal. Also, the price of an American lookback put is always 
greater than the European one. 

However, table 4 shows that the upper bounds on the lookback put price derived with Snell 
envelopes techniques in Conze and Viswanathan (1991), are quite loose. This is especially 
true for short maturities, for which the price is overestimated by a factor of more than 100% 
on average. 

8.2 Average-rate option 

Results for fixed-strike and floating-strike average-rate options are presented in tables 6 and 5 
respectively. Because the average variable needs a dense quantization (p = 0.1), the required 
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Lookback call option 
S 0 = 100, r = 10% 





European 


a 


T 


analytic 


step 
(day) 


FSG 

C 


(C a - C)/C 












2 


8.250 


7.71% 








1 


8.480 


5.13% 


20 


3 


8.939 


0.5 


8.611 


3.67% 








0.1 


8.790 


1.66% 








0.05 


8.833 


1.12% 



Table 3: Measured convergence of the FSG method for a lookback option. 



Lookback put option 
S 0 = 100, r = 10% 





American 


a 


T 


FSG Snell u. b. 
P P 


error 

\(P-P)/P\ 






10 


3 


3.087 5.260 


70% 


6 


4.209 8.653 


105 % 


12 


5.630 14.920 


165 % 






20 


3 


6.845 9.121 


33 % 


6 


9.606 14.288 


49 % 


12 


13.307 23.135 


74% 






40 


3 


14.811 17.294 


17% 


6 


21.362 26.614 


25 % 


12 


30.661 42.067 


37% 



Table 4: Comparison of American lookback FSG put prices and corresponding Snell envelope 
upper bounds. 
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Floating-strike average-rate call option 
So = 100, r = 10% 





European 


American 


a 


T 


C 


10 


3 


1.852 


2.142 


6 


3.092 


3.404 


12 


5.382 


5.670 






20 


3 


2.943 


3.689 


6 


4.555 


5.514 


12 


7.259 


8.463 






40 


3 


5.195 


6.817 


6 


7.673 


9.860 


12 


11.509 


14.446 



Floating- strike average-rate put option 
So = 100, r = 10% 





European 


American 


a 


T 


P 


10 


3 


0.632 


1.194 


6 


0.671 


1.494 


12 


0.614 


1.799 






20 


3 


1.724 


2.773 


6 


2.135 


3.684 


12 


2.498 


4.812 






40 


3 


3.981 


5.996 


6 


5.260 


8.223 


12 


6.785 


11.229 



Table 5: Results of the FSG method for a floating-strike average-rate option. Time Steps were 
set to 1, 2, and 3 days for 3, 6, and 12 months maturities respectively. Precision for FSG prices 
is 0.1%. 



memory is l/p = 10 times as much as in the lookback case. Maintaining a 1 day time step for 
all maturities would thus require too much memory on casual machines. Time Steps were set to 
1,2, and 3 days for 3, 6, and 12 months maturities respectively, so as to get about N = 100 time 
steps in each case. This is still quite reasonable. Time and space complexities are 0(N 3 /p) 
and 0(N 2 / p) respectively. The computation time for one call/put is about 15 second on a DEC 
alpha PC, and the required memory 4 Megabytes. 

There is no analytic solution to the average-rate option pricing problem in general. However, 
as seen in sec. 3.3, the zero-strike case leads to the closed-form formula of eq. (3.6). Table 7 
compares the analytic prices, and those computed with the FSG method in this case. The preci- 
sion is 0.1% at least. Because the PDE (3.5) governing the average-rate price is homogeneous 
in the spot price So and the strike K, the convergence of the FSG method does not depend on 
the value of K. Therefore, given the quantization with p = 0.1 of A, we can state an overall 
precision of at least 0.1% in our results. 

Average-rate options are much cheaper than ordinary and lookback options. Also, since the 
average is not a strictly decreasing (increasing) time process as is the case for the minimum 
(maximum), American average -rate calls are worth exercising. 

A possible method for approximating European call and put prices C ar ith and P ar ith of a fixed- 
strike average-rate option consists in replacing the arithmetic time-average with the geometric 
time-average, in which case closed-form formulas can be derived for prices C geo and P geo 
(Kemna and Vorst (1990)). Since the geometric average is always less than the arithmetic 
average, the following relationships holds: C geo < C ar ith, and P geo > P ar ith- 
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Fixed-strike average-rate call option 
So = 100, r = 10% 





European 


American 


a 


T 


K 


C 


10 


3 


95 

100 

105 


6.132 
1.869 
0.151 


6.546 
1.967 
0.152 




6 


95 

100 

105 


7.248 
3.100 
0.727 


7.632 
3.212 
0.735 




12 


95 

100 

105 


9.313 
5.279 
2.313 


9.616 
5.394 
2.336 






20 


3 


95 

100 

105 


6.500 
2.960 
0.966 


7.371 
3.219 
1.001 




6 


95 

100 

105 


7.793 
4.548 
2.241 


8.805 
4.893 
2.337 




12 


95 

100 

105 


10.336 
7.079 
4.539 


11.218 
7.521 
4.729 






40 


3 


95 

100 

105 


8.151 
5.218 
3.106 


9.447 
5.826 
3.347 




6 


95 

100 

105 


10.425 
7.650 
5.444 


10.927 
8.519 
5.913 




12 


95 

100 

105 


13.825 
11.213 
8.989 


15.649 
12.439 
9.790 



Fixed-strike average-rate put option 
So = 100, r = 10% 





European 


American 


a 


T 


K 


P 






95 


0.013 


0.013 




3 


100 


0.626 


0.832 






105 


3.785 


5.337 










95 


0.046 


0.051 


10 


6 


100 


0.655 


0.978 






105 


3.039 


5.287 










95 


0.084 


0.104 




12 


100 


0.577 


1.079 






105 


2.137 


5.230 










95 


0.379 


0.407 




3 


100 


1.716 


2.066 






105 


4.598 


6.108 










95 


0.731 


0.820 


20 


6 


100 


2.102 


2.629 






105 


4.552 


6.338 










95 


1.099 


1.318 




12 


100 


2.369 


3.181 






105 


4.356 


6.596 










95 


2.025 


2.223 




3 


100 


3.970 


4.581 






105 


6.735 


8.168 










95 


3.215 


3.610 


40 


6 


100 


5.197 


6.078 






105 


7.748 


9.438 










95 


4.550 


5.263 




12 


100 


6.465 


7.761 






105 


8.767 


10.927 



Table 6: Results of the FSG method for a fixed-strike average-rate option. Time Steps were 
set to 1, 2, and 3 days for 3, 6, and 12 months maturities respectively. Precision for all prices 
is 0.1%. 
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Zero-strike average-rate call option 
S 0 = 100, r = 10% 





European 


a 


T 


analytic FSG 

C a C 


\(C a -C)/C\ 






10 


3 


98.763 98.780 


0.02% 


6 


97.547 97.579 


0.03% 


12 


95.175 95.232 


0.05% 






20 


3 


98.763 98.781 


0.02% 


6 


97.547 97.580 


0.03% 


12 


95.175 95.239 


0.07% 






40 


3 


98.763 98.785 


0.02% 


6 


97.547 97.587 


0.04% 


12 


95.175 95.277 


0.10% 



Table 7: Accuracy of the FSG method measured against the zero-strike average-rate option. 
Time Steps were set to 1, 2, and 3 days for 3, 6, and 12 months maturities respectively. 

Table 8 shows that the geometric average approximation is good for low terms and low 
volatilities. However, the error rapidly grows with the option term and volatility, and stops 
being acceptable (up to 10%) for one year and longer standing options. 

9 Conclusion 

In this paper, we have analyzed the problem of pricing path-dependent contingent claims. 
We have shown that these problems lead to solving degenerate diffusion PDE in the space 
augmented with the path-dependent variables. We have established necessary and sufficient 
conditions under which these degenerate PDE are holonomic, i.e. can be reduced to lower 
dimensional non-degenerate PDE. We have applied these results to popular types of path- 
dependent options. In particular, we have shown that the average-rate option pricing problem 
is non-holonomic. 

Then, we have described a new numerical technique called the Forward Shooting Grid method 
(FSG) for pricing both European and American non-holonomic contingent claims, and have 
tested it on lookback and average-rate options. It is straightforward to use the FSG method in 
order to price other popular path-dependent options, such as capped or barrier options. 

The FSG method proves to be as accurate as Monte Carlo simulation, with faster execution 
time. It is also the first method capable of dealing with the early exercise condition of American 
path-dependent options, showing that the Snell envelope upper bounds obtained on American 
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Geometric average-rate call option 
S 0 = 100, r = 10% 





European 


o 


T 


K 


analytic 


error w.r.t. 
FSG arith 


10 


3 


95 

100 

105 


6.093 
1.830 
0.141 


0.63 % 
2.08 % 
6.62 % 




6 


95 

100 

105 


7.165 
3.022 
0.683 


1.14% 
2.52 % 
6.05 % 




12 


95 

100 

105 


9.160 
5.135 
2.196 


1.64% 
2.73 % 
5.06% 






20 


3 


95 

100 

105 


6.402 
2.872 
0.908 


1.51 % 
2.97 % 
6.00 % 




6 


95 

100 

105 


7.742 
4.374 
2.102 


0.65 % 
3.83 % 
6.20 % 




12 


95 

100 

105 


9.979 
6.751 
4.251 


3.45 % 
4.63 % 
6.34 % 






40 


3 


95 

100 

105 


7.862 
4.959 
2.889 


3.54% 
4.96 % 
6.99 % 




6 


95 

100 

105 


9.867 
7.136 
4.983 


5.35 % 
6.72 % 
8.47 % 




12 


95 

100 

105 


12.743 
10.200 
8.048 


7.83 % 
9.03 % 
10.47 % 



Geometric average-rate put option 
So = 100, r = 10% 





European 


(7 


r 


K 


analytic 


error w.r.t. 
FSG arith 


10 


3 


95 

100 

105 


0.013 
0.627 
3.815 


-0.00 % 
-0.16% 
-0.80 % 




6 


95 

100 

105 


0.047 
0.660 
3.079 


-2.17% 
-0.76 % 
-1.32% 




12 


95 

100 

105 


0.086 
0.587 
2.173 


-2.38 % 
-1.73 % 
-1.68% 






20 


3 


95 

100 

105 


0.384 
1.730 
4.643 


-1.32% 
-0.82 % 
-0.98 % 




6 


95 

100 

105 


0.746 
2.134 
4.619 


-2.05 % 
-1.52% 
-1.47% 




12 


95 

100 

105 


1.142 
2.440 
4.465 


-3.91 % 
-3.00 % 
-2.50 % 






40 


3 


95 

100 

105 


2.089 
4.064 
6.870 


-3.16% 
-2.37 % 
-2.00 % 




6 


95 

100 

105 


3.355 
5.380 
7.985 


-4.35 % 
-3.52 % 
-3.06 % 




12 


95 

100 

105 


4.847 
6.829 
9.203 


-6.53 % 
-5.63 % 
-4.97 % 



Table 8: Comparison between the European prices of a fixed-strike geometric average-rate call 
and a fixed-striked (arithmetic) average-rate call. FSG prices, which are not reported here, 
were taken from table 6 in order to compute the relative error. 



Research Report No. 37 



January 1994 



36 



Jerome Barraquand and Thierry Pudet 



lookback prices are quite overestimated. Our numerical experiments have also shown that the 
usual geometric average approximation of arithmetic average-rate options is inaccurate. 

The FSG method is a general purpose solution technique for arbitrary multidimensional 
advection-diffusion equations. However, since these problems require a memory space ex- 
ponential in the number of variables, the FSG method, like any finite difference or lattice-based 
method, can only be used on problems with few variables. 

Unlike Finite Difference Methods, the FSG method is unconditionally convergent, even when 
the diffusion term is degenerate. Such degeneracy arises in several other important pricing 
problems. In particular, we plan to investigate the application of the FSG method to the pricing 
of path-dependent interest rate contingent claims such as Mortgage-Backed securities. 
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